<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Minimize stopband ripple of a linear phase lowpass FIR filter</title>
<link rel="canonical" href="/Users/mcgrant/Projects/CVX/examples/filter_design/html/fir_lin_phase_lowpass_min_ripple.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Minimize stopband ripple of a linear phase lowpass FIR filter</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs a linear phase FIR lowpass filter such that it:</span>
<span class="comment">% - minimizes the maximum passband ripple</span>
<span class="comment">% - has a constraint on the maximum stopband attenuation</span>
<span class="comment">%</span>
<span class="comment">% This is a convex problem.</span>
<span class="comment">%</span>
<span class="comment">%   minimize   delta</span>
<span class="comment">%       s.t.   1/delta &lt;= H(w) &lt;= delta     for w in the passband</span>
<span class="comment">%              |H(w)| &lt;= atten_level        for w in the stopband</span>
<span class="comment">%</span>
<span class="comment">% where H is the frequency response function and variables are</span>
<span class="comment">% delta and h (the filter impulse response).</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% user's filter specifications</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% filter order is 2n+1 (symmetric around the half-point)</span>
n = 10;

wpass = 0.12*pi;        <span class="comment">% passband cutoff freq (in radians)</span>
wstop = 0.24*pi;        <span class="comment">% stopband start freq (in radians)</span>
atten_level = -30;      <span class="comment">% stopband attenuation level in dB</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% create optimization parameters</span>
<span class="comment">%********************************************************************</span>
N = 30*n+1;                            <span class="comment">% freq samples (rule-of-thumb)</span>
w = linspace(0,pi,N);
A = [ones(N,1) 2*cos(kron(w',[1:n]))]; <span class="comment">% matrix of cosines</span>

<span class="comment">% passband 0 &lt;= w &lt;= w_pass</span>
ind = find((0 &lt;= w) &amp; (w &lt;= wpass));   <span class="comment">% passband</span>
Ap  = A(ind,:);

<span class="comment">% transition band is not constrained (w_pass &lt;= w &lt;= w_stop)</span>

<span class="comment">% stopband (w_stop &lt;= w)</span>
ind = find((wstop &lt;= w) &amp; (w &lt;= pi));  <span class="comment">% stopband</span>
Us  = 10^(atten_level/20)*ones(length(ind),1);
As  = A(ind,:);

<span class="comment">%********************************************************************</span>
<span class="comment">% optimization</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% formulate and solve the linear-phase lowpass filter design</span>
cvx_begin
  variable <span class="string">delta</span>
  variable <span class="string">h(n+1,1)</span>;

  minimize( delta )
  subject <span class="string">to</span>
    <span class="comment">% passband bounds</span>
    Ap*h &lt;= delta;
    inv_pos(Ap*h) &lt;= delta;

    <span class="comment">% stopband bounds</span>
    abs( As*h ) &lt;= Us;
cvx_end

<span class="comment">% check if problem was successfully solved</span>
disp([<span class="string">'Problem is '</span> cvx_status])
<span class="keyword">if</span> ~strfind(cvx_status,<span class="string">'Solved'</span>)
  <span class="keyword">return</span>
<span class="keyword">else</span>
  <span class="comment">% construct the full impulse response</span>
  h = [flipud(h(2:end)); h];
  fprintf(1,<span class="string">'The optimal minimum passband ripple is %4.3f dB.\n\n'</span>,<span class="keyword">...</span>
            20*log10(delta));
<span class="keyword">end</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% plots</span>
<span class="comment">%********************************************************************</span>
figure(1)
<span class="comment">% FIR impulse response</span>
plot([0:2*n],h',<span class="string">'o'</span>,[0:2*n],h',<span class="string">'b:'</span>)
xlabel(<span class="string">'t'</span>), ylabel(<span class="string">'h(t)'</span>)

figure(2)
<span class="comment">% frequency response</span>
H = exp(-j*kron(w',[0:2*n]))*h;
<span class="comment">% magnitude</span>
subplot(2,1,1)
plot(w,20*log10(abs(H)),[wstop pi],[atten_level atten_level],<span class="string">'r--'</span>);
axis([0,pi,-40,10])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'mag H(w) in dB'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling SDPT3 4.0: 872 variables, 278 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 278
 dim. of sdp    var  = 74,   num. of sdp  blk  = 37
 dim. of socp   var  = 458,   num. of socp blk  = 229
 dim. of linear var  = 303
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|7.3e+02|4.9e+01|1.1e+05| 2.048238e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|0.921|0.699|5.8e+01|1.5e+01|2.5e+04| 5.008791e+01 -2.175738e+01| 0:0:00| chol  1  1 
 2|0.732|0.870|1.6e+01|1.9e+00|4.3e+03| 6.281441e+01 -4.143755e+01| 0:0:00| chol  1  1 
 3|0.995|0.863|7.8e-02|2.7e-01|5.4e+02| 6.205041e+01 -4.108188e+01| 0:0:00| chol  1  1 
 4|1.000|0.942|5.0e-04|1.6e-02|7.0e+01| 4.728715e+01 -3.430874e+00| 0:0:00| chol  1  1 
 5|0.368|1.000|3.1e-04|1.1e-04|4.3e+01| 3.624702e+01 -6.510807e+00| 0:0:00| chol  1  1 
 6|1.000|0.921|1.5e-05|7.2e-05|6.5e+00| 3.973191e+00 -2.473693e+00| 0:0:00| chol  1  1 
 7|0.916|0.900|1.4e-06|1.0e-05|2.8e+00| 1.071821e+00 -1.754710e+00| 0:0:00| chol  1  1 
 8|0.820|0.946|2.6e-07|8.4e-07|1.0e+00|-1.829241e-01 -1.179210e+00| 0:0:00| chol  1  1 
 9|1.000|1.000|1.6e-09|5.4e-08|3.8e-01|-7.085424e-01 -1.089012e+00| 0:0:00| chol  1  1 
10|0.718|0.879|4.5e-10|6.9e-09|1.4e-01|-9.217470e-01 -1.059581e+00| 0:0:00| chol  1  1 
11|0.958|0.914|1.9e-11|7.0e-10|1.9e-02|-1.034838e+00 -1.054108e+00| 0:0:00| chol  1  1 
12|0.819|0.737|3.5e-12|1.9e-10|7.8e-03|-1.044891e+00 -1.052643e+00| 0:0:00| chol  1  1 
13|0.795|0.928|6.9e-13|1.5e-11|3.0e-03|-1.048933e+00 -1.051939e+00| 0:0:00| chol  1  1 
14|0.662|0.907|2.0e-13|2.4e-12|1.5e-03|-1.050281e+00 -1.051737e+00| 0:0:00| chol  1  1 
15|0.711|0.932|2.7e-12|1.2e-12|6.4e-04|-1.050991e+00 -1.051632e+00| 0:0:00| chol  1  1 
16|0.787|0.942|5.0e-12|1.1e-12|2.3e-04|-1.051364e+00 -1.051591e+00| 0:0:00| chol  1  1 
17|0.929|0.860|1.7e-12|1.1e-12|3.9e-05|-1.051542e+00 -1.051581e+00| 0:0:00| chol  1  1 
18|0.910|0.810|1.4e-12|1.2e-12|5.9e-06|-1.051573e+00 -1.051579e+00| 0:0:00| chol  1  1 
19|0.897|0.865|1.6e-12|1.2e-12|8.9e-07|-1.051577e+00 -1.051578e+00| 0:0:01| chol  1  1 
20|0.945|0.969|2.1e-12|1.0e-12|7.6e-08|-1.051578e+00 -1.051578e+00| 0:0:01| chol  1  1 
21|0.725|1.000|7.0e-13|1.0e-12|3.9e-08|-1.051578e+00 -1.051578e+00| 0:0:01|
  stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 21
 primal objective value = -1.05157797e+00
 dual   objective value = -1.05157801e+00
 gap := trace(XZ)       = 3.93e-08
 relative gap           = 1.27e-08
 actual relative gap    = 1.27e-08
 rel. primal infeas (scaled problem)   = 7.01e-13
 rel. dual     "        "       "      = 1.00e-12
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.2e+00, 1.1e+00, 1.2e+01
 norm(A), norm(b), norm(C) = 8.6e+01, 2.0e+00, 9.6e+00
 Total CPU time (secs)  = 0.57  
 CPU time per iteration = 0.03  
 termination code       =  0
 DIMACS: 7.0e-13  0.0e+00  4.8e-12  0.0e+00  1.3e-08  1.3e-08
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.05158
 
Problem is Solved
The optimal minimum passband ripple is 0.437 dB.

</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_lin_phase_lowpass_min_ripple__01.png" alt=""> <img src="fir_lin_phase_lowpass_min_ripple__02.png" alt=""> 
</div>
</div>
</body>
</html>